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We solve the dynamical GR equations for the spherically symmetric evolution 
of compact stars in the vicinity of the maximum mass, for which instability sets in 
according to linear perturbation theory. The calculations are done with the analytical 
Zeldovich-like EOS P = a{p — po) and with the TMl parametrisation of the RMF 
model. The initial configurations for the dynamical calculations are represented by 
spherical stars with equilibrium density profile, which are perturbed by either (i) 
an artificially added inward velocity field proportional to the radial coordinate, or 
(ii) a rarefaction corresponding to a static and expanded star. These configurations 
are evolved using a one-dimensional GR hydro code for ideal and barotropic fluids. 
Depending on the initial conditions we obtain either stable oscillations or the collapse 
to a black hole. The minimal amplitude of the perturbation, needed to trigger 
gravitational collapse is evaluated. The approximate independence of this energy 
on the type of perturbation is pointed out. At the threshold we find type I critical 
behaviour for all stellar models considered and discuss the dependence of the time 
scaling exponent on the baryon mass and EOS. 
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I. INTRODUCTION 


The research on stellar oscillations within the framework of GR goes back to the early 
works of Chandrasekhar and others from about half a century ago The most simple 

class of oscillations are radial oscillations, which can be classihed by a single index denoting 
the number of nodes within the star. Within the sequence of equilibrium configurations 
for a given equation of state (EOS), there is one for which the squared frequency of the 
fundamental mode crosses zero. It is usually identified as the critical star which corresponds 
to the configuration of maximum mass and maximum baryon number. In the mass-radius 
diagramm it separates the branch of perturbatively unstable stars from the stable branch. 

The dynamics of gravitational collapse in GR was studied for the hrst time in j^, 
see also references therein. The authors considered a polytropic equation of state P = KrP 
and a large number of different choices for 7, the total mass and the radius of the initial 
conhgurations, which were assumed to be at rest. The perturbations were imposed by 
streching the matter from its equilibrium radius to some larger radius and thereby obtaining 
an inner pressure depletion. Specihcally, the authors asked the question, whether black 
holes could be formed out of stellar cores with a mass below the maximum mass for a given 
EOS. The conclusion was, that stars below the maximum mass are dynamically stable even 
for large perturbations, i.e. they perform oscillations without evolving to a black hole. 


Similar studies with a more realistic microphysical description of the EOS were per¬ 
formed in 6|, l7j . The initial configurations were neutron stars with an equilibrium density 
prohle, to which an inward linear velocity held was added by hand. In contrast to it was 
found that linearly stable neutron stars can collapse to a black hole, if the perturbation is 
large enough and the mass of the star exceeds a lower threshold depending on the EOS. 

Non-linear effects in the evolution of spherical stars have been studied in 8|, l9(]. A 
detailed investigation of the nonlinear dynamics and its phenomenology was obtained 
in 1^. It was found that mode pairs with an integer frequency ratio can efficiently 
interact and exchange energy. Furthermore, it was shown that linearly unstable stars can 
perform stable oscillations over a long period of time without collapsing to a black hole. 
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The evolution of linearly super-critical stars in GR was also investigated in [n|. They 
considered polytropes with a constant baryon mass and showed that the lifetime of the 
metastable solution along the evolutionary path exhibits the scaling relationship 


r = —a In \p — p*\ + const 


( 1 ) 


, which dehnes the critical phenomenon of type I 12|. Here, a denotes the time scal¬ 
ing exponent, which controls the sharpness of the threshold behaviour near black hole 
formation. The metastable escape time or lifetime of the intermediate state r measures 
how long the solution with parameter p stays close to the critical solution with parameter p*. 


A thorough investigation of the critical collapse dynamics of low-mass neutron stars 
has been presented in tl3|, where the scaling behaviour of black hole formation was 

npn 

studied. This letter aims to extend the previous work in [5H7|, 113| and to investigate radial 
oscillations of large amplitudes in the vicinity of the maximum-mass conhguration, taking 
into account non-linear effects. Specihcally, we want to address the question, how strong 
the perturbation should be, in order to traverse the barrier and cause a collapse to a black 
hole (see £g. 7). In contrast to previous work, we will identify conhgurations by their 
baryon number and therefore avoid the decomposition into background and perturbation, 
which is arbitrary to some extent. 


It is clear that the energy barrier, which separates stable neutron stars from black 
holes, must change smoothly in the vicinity of the maximal mass. In this letter we report 
on a systematic survey of this energy wall of massive compact stars and its dependence 
on EOS stiffness and baryon number, which is lacking in previously published work. The 
envisioned collapse may be triggered by mass accretion in binary systems or in a supernova 
explosion, where the expelled layers fail to reach the escape velocity and fall back onto the 
proto-neutron star. 
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II. DYNAMICAL EQUATIONS 


In this section we consider non-linear radial oscillations of compact stars. We keep the 
assnmption of spherical symmetry, bnt allow for perturbations of arbitrary amplitude. We 
assume that there are no gradients in the specihc entropy (homentropic fluid). In this case 
the equation for baryon number conservation is not linearly independent from the equation 
for energy-momentum conservation. The equations governing the dynamics are Einstein’s 
equation and the energy-momentum conservation equation.^ 

( 2 ) 

= 0 (3) 

We adopt the Schwarzschild-like parametrisation of the metric 

ds^ = ^ e2A(t,r)^^2 ^ (4) 


with time-dependent metric functions <I> and A. The normalisation of the 4-velocity yields 

^2A-2<I>„,i2 2 $ 




u -\- e 


( 5 ) 


a relationship between and (henceforth denoted by u). It should be stressed that our 
choice of the metric involves a coordinate singularity at the Schwarzschild radius and thus 
we cannot follow the motion of matter trough the event horizon. We introduce the rescaled 
radial coordinate rj = r/R{t), where R{t) is the time-dependent stellar radius. Then eq. ([2]) 
leads to the two constraint equations for the metric functions 

p2A _ 1 

= AttR^tj ip + P) -I- AnR'^riPe^^ H- (6) 

2t] 

1 - 

A,„ = AttR^t] ip + P) -|- AirR^ppe^^ -\ - (7) 

2p 

where denotes a partial derivative, p the total energy density (including rest mass) and 
P the pressure. Adopting the energy-momentum tensor for an ideal fluid, eq. (|3]) provides 
a hyperbolic set of two evolution equations for the pressure and fluid 4-velocity 


-P,t + CiiP,ri +Ci2U,r) +C13 — 0 


( 8 ) 


^ Here and below we use geometrical units, G = c = km = 1. 
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U,t + C2lP,r] +CiiU,ri +C23 — 0 


(9) 


Here we have introduced the notations: 

(i -1) 


Cll = 




Rf] 

'r 


C12 = 


p + P 

Ru% 


$ 0 
e u u 


Cl3 = 


C21 = 


C 23 — 


2 r (1 + i _ 

dp 2 »I»- 2 A „,0 

^ 

R{p + P)i 
,$-2A.,0 r Qp 


Se'*’ (p + P) + (p + P) ( — 1 — S'liPr'^ + 4-u^) 


u 


2ri 


and 


^ (e'*’ + (-1 - SvrPr^) j + 4e 




2A+4'^2 


( 10 ) 

( 11 ) 

( 12 ) 

( 13 ) 

( 14 ) 


( 16 ) 


The set of eqs. for the hve unknown variables P, p, <h, A is given by eqs. l|6|) l|7|) (|8|1 (pjl 
together with the equation of state. The enclosed baryon number Nb obeys the differential 
equation 

Nb ,r, = 47rP^pne^+^M°, (16) 

where n denotes the baryon number density. The evolution of the star radius is obtained 
from the equation 

u{R) 


R = 


u\R)' 


(17) 


Equations ([8]) and ([9]) are equivalent to eqs. (5.85) and (5.86) in js], as can be demonstrated 
by transforming the radial coordinate from p to r. It is easy to show that for the stationary 
limit eq. ([8]) is trivially fullhlled, whereas eq. (|9]) leads to the well known TOV equation of 

equation for radial modes can be obtained by 


relativistic stellar structure. The eigenmode 

nn 

perturbing eq. (|9]) around equilibrium [Jj, [2|, and imposing harmonic time-dependence. 
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III. EQUATION OF STATE 


A. Zeldovich EOS 


In our calculations we have used two equations of state. The hrst one is the Generalised 
Zeldovich (GZ) EOS P = a{p — po), which for a = 1 yields the causal-limit EOS proposed 
by Ya. Zeldovich Q and which for a = 1/3 coincides with the MIT bag model jl^. This 
simple parametrisation is able to qualitatively reproduce the bulk properties of neutron 
stars. While a controls the stiffness and the sound velocity, the parameter po introduces a 
dimensional scale in the EOS. Different choices for a and po allow us to consider families 
with different maximum mass or radius. Interestingly, the compactness (= 2 M/R) of the 
maximum mass or maximum radius conhguration does not depend on po. In Eg. 1 we plot 
the mass of the maximum-mass conhguration for the parameter space given by a and po. 
The high mass stars in the bottom-right corner correspond to low central densities, while the 
low mass stars in the top-left corner correspond to large central densities. In Eg. 2 we show 
the radii of the same stellar models as in Eg. 1. The baryon number density n is related to 
the energy density p by n/riQ = [a (p -|- jqgj-g jg ^ scaling variable, which 

needs to be found from a speciEc matter model. We note that the baryon number ratio of 
two stellar models does not depend on the choice of uq. This equation solves the well known 
identity for the baryon chemical potential pb of cold neutron star matter 


dn n 

dp p-|-P 


(18) 


B. RMF EOS 


For a realistic description of neutron stars, the relativistic mean held models are often 
employed 19(]. They provide an efiective and relativistically covariant description of dense 
baryonic systems. The parameters are Etted to the known properties of nuclear matter 
at saturation. As the second EOS we use the TMl parameter set 20|, which leads to a 
maximum mass of about 2.2Mq. 
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FIG. 1: Maximum mass M of compact stars with the generalised Zeldovich EOS P = a{p — po) 
as function of a and po. The lines indicate configurations of constant M/Msun- The choices of a, 
Po for the simulations are indicated by solid squares. 



FIG. 2: Radius R of the maximum mass configuration of compact stars with the generalised 
Zeldovich EOS P = a {p — po) as function of a and po- The lines indicate configurations of constant 
R/km. 











































IV. RESULTS 


We consider 3 different classes of star configurations based on the above-mentioned EOS, 
namely 

(i) the GZ model with a = 1, po = 7.5367 

(ii) the GZ model with a = 1/3, po = 2.7782 

(iii) the TMl EOS based on the relativistic mean held model. 

The parameters in the hrst two cases are chosen to yield a maximum mass of precisely 2.1Mq 
(see hg. 1), whereas the TMl EOS yields a maximum mass of approximately 2.2Mq. For 
each EOS we select 10 models based on the baryon number of the maximum mass conhgu- 
ration NB,max, namely conhgurations with a baryon number Nb = {1 — k 10“^) NB,max with 
k = 1,2..10. 


A. Initialisation 

We consider two different types of initial perturbations preserving the baryon number, 
namely the velocity held induced perturbation and the rarefaction induced perturbation. 
Velocity field induced perturbation-. Once Nb is specihed, we choose the central energy 
density, which leads to a stable equilibrium solution with baryon number Nb- We then 
impose an inward velocity held uirj) = —p Umax, where Umax is the surface velocity. To 
preserve Nb according to eq. (ITOll . the density and pressure prohles of the TOV solutions 
have to be adapted. We therefore select a radius independent x and redehne the pressure 
according to Pnewiji) = Poidiji) (1 + ^:). Since the EOS is barotropic, this procedure also 
determines the energy density. This rescaling allows to preserve Nb, while considering 
perturbations of diherent strength. Note that this procedure does not ahect the stellar 
radius. 

Rarefaction induced perturbation-. Once Nb is selected, we again choose the central energy 
density, which leads to a stable equilibrium solution with baryon number Nb- The stellar 
radius is then artihcially increased, while p {p) and P {p) are held constant. In order to 
preserve Nb according to eq. (USD, we apply the same rescaling of density and pressure 
as described above. In contrast to above, here the initial stellar radius is not preserved 
and the perturbed star is initially at rest. Despite the names we attach to the two distinct 
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time (milli seconds) 


FIG. 3: The time evolution of the stellar radius for the TMl EOS with a baryon number 4% 
below the maximum value. 4 trajectories corresponding to a different strength of the velocity field 
induced perturbation are shown. The velocity profile belonging to the red solid line is shown in 
fig. 5. 




FIG. 4; The metastable escape time (left) and its relative error with respect to the linear ht (right) 
as function of initial data. The legend indicates the specific stellar models, namely 1: (iii,l,v), 2: 
(iii,I0,v), 3: (iii,10,r), 4: (i,10,v), 5: (iii,4,r), 6: (ii,4,r). 
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FIG. 5: The radial velocity u as a function of time and radius for the slightly supercritical solution 
shown in fig. 3. The green contours denote zero velocity points and separate infalling matter from 
expanding matter. 
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FIG. 6: The radial velocity u as a function of time and radius for a slightly supercritical solution 
belonging to the case (ii,4,v). The green contours denote static fluid elements. 
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kinds of initial pertnrbations, both involve a rarefaction, which in the hrst case however is 
considerably weaker. Henceforth we use tripels (EOS type {i,ii,iii}, baryon dehcit {1,2..10}, 
perturbation type {v,r}) to label each of the 60 cases considered in the present study. 

Figure 3 shows the evolution of the stellar radius for different magnitudes of the ve¬ 
locity induced perturbations. If the excitation energy is subcritical, the star performs 
oscillations of the fundamental mode, together with overtones. As the energy is shifted 
closer to criticality, the frequency of the fundamental oscillation decreases. We also notice 
that the shape of the oscillations deviates more and more from the harmonic curve. 
The closer the star is to the critical regime, the longer are delays at the compression 
phase. In other words, the closeness to criticality manifests itself in a long frustrated 
phase, which is characterised by very small amplitude oscillations and can extend over 
several milli seconds. If the excitation energy is supercritical, the star evolves to a black hole. 

For each case the bisection method is used to tune the strength of the initial pertur¬ 
bation towards criticality at the threshold between secular oscillations and gravitational 
collapse. Figure 4 shows the escape time r and its deviation from a linear £t as function of 
the closeness of initial data to criticality. The parameter p has been identihed with u^ax for 
velocity induced perturbations and with the initial excess in radius AR otherwise. Here, r 
has been identihed as the length of the time interval during which the stellar radius deviates 
by less than 2% from the solution obtained at the last bisection step.^ For any value of p* 
a unique time scaling relation ensues. To determine the critical p* in the sense of eq. ([I]) 
we minimise the error of the linear £t of the ensuing relation with respect to a free p*.^ 
The energy associated with this perturbation is denoted as critical excitation energy or 
short critical energy. It is the energy difference between the critical solution and the stable 
equilibrium solution of the same baryon number and can be calculated from eq. O- 

The result of our analysis is shown in £g. 7, where we plot the critical excitation 

energy as function of the baryon number in units of the maximum baryon number of the 

^ Only subcritical evolutions have been used for the scaling relationship. The inclusion of supercritical ones 
does not change the picture. 

^ In practice this yields a value very close to the last bisection step. 
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selected EOS. The critical excitation energy decreases with increasing baryon number 
and the curves converge as the distance to the maximum mass conhguration diminishes. 
Furthermore, regardless of the baryon number considered, the soft GZ case (case ii) exhibits 
a smaller critical energy than the TMl case (case hi), which itself lies below the stiff GZ case 
(case i). The dependence on the type of perturbation is remarkably weak. For the case of the 
velocity induced perturbation we hnd an invariance of the critical energy with respect to the 
sign of the velocity, i.e. the result is exactly the same for initially ingoing and outgoing flow. 

Next, we study the behaviour of the time scaling exponent a as function of the 
baryon number. We recall from eq. ([I]) that a is dehned as the negative slope of the curve, 
which relates the lifetime of the intermediate state with In \p — p*\ for any one-dimensional 
parametetrisation p of initial data. Within the parameter space of initial data, a controls 
the width of the neighbourhood around criticality, for which threshold behaviour is 
signihcant. In practice, we have obtained a as the negative slope of the linear fit for the 
subset of subcritical evolutions. 


The result of our analysis is shown in hg.8. The scaling exponent steadily increases 
with baryon number and gains about 70% within the band of rest masses considered. 
There is no noticeable dependence on the perturbation type. Remarkably, we hnd identical 
behaviour for the EOS cases (ii) and (iii), i.e. strange stars with the simple MIT bag model 
and hadronic stars with the TMl EOS exhibit exactly the same time scaling exponent. 
In contrast, for case (i) we hnd a reduction of a by about 30%. Our reported values 
consistently exceed those reported in [2l| for polytropic ’’Gaussian packet systems” by a 
factor of 1.2 to 2.5. We hnd the same dependence on baryon number, even though our 
result seems to support a convex dependence and we cannot conhrm their battening at high 
baryon masses. 


V. CONCLUSION 

In this work we have studied the dynamical evolution of nonlinear perturbations on 
spherical conhgurations with a baryon mass above 90% of the maximal value of the respec¬ 
tive EOS. Numerical simulations in full GR were performed with no other restrictions, than 
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GO a= 1, velocity field 
AAa=1 , rarefaction 
WTM1, velocity field 
-H- TM1, rarefaction 
xx a=1/3, velocity field 
□□a= 1/3, rarefaction . 
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FIG. 7: The critical excitation energy as function of the baryon number. The unit on the x-axis 
is the maximum baryon number of the considered EOS, which is arbitrary for case (i) (see free 
scaling variable above), 3.4773 10®^ for case (ii) and 3.0704 10^'^ for case (hi). 
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FIG. 8; The time scaling exponent (see eq. ([T])) as function of the baryon number for the several 
considered EOS and perturbation types. 
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spherical symmetry and ideal hydrodynamics. The calcnlations were done for 3 different 
EOS and 10 different valnes of the stellar rest mass. A bisection strategy was invoked 
to determine the minimal energy for each EOS and baryon mass, which leads to black 
hole formation. This energy tnrns ont to be more than an order of magnitnde below the 
mass difference to the maximnm mass conhgnration of the considered EOS. The critical 
energy increases with decreasing baryon nnmber, i.e. the stars on the stable branch of the 
mass-radins diagram are stable to a different extent. Fnrthermore its dependence on the 
type of pertnrbation is weak, bnt we have fonnd a noticeable dependence on the EOS. 
Compact stars with stiff eqnations of state tend to be more resistent to collapse than their 
soft connterparts. In particnlar, strange qnark stars shonld collapse easier to black holes 
than their hadronic connterparts. 

Fnrthermore, we have analysed the threshold behavionr close to black hole formation 
and fonnd excellent agreement with the scaling relationship ([T]) for all the models consid¬ 
ered. We therefore conjectnre, that type 1 critical behavionr is a generic featnre of high-mass 
compact stars. Remarkably, we have fonnd nniversal scaling behavionr for hadronic and 
strange qnark stars, independent of the pertnrbation type. On the other hand, the weak 
dependence of the excitation energy on the specihc type of initial pertnrbation snggests 
that the critical solntion is not nniversal with respect to different families of initial data 
and hence the initialisation of the pertnrbations in GR-critical collapse scenarios deserves 
particnlar attention. Whether this is dne to an inconsistent treatment of shocks and 
thermodynamics in onr simnlations of the lower mass stars or a non-nniversality of the 
critical solntion, we wish to explore in fntnre work. 
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